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Abstract 

We determine the spectrum of B s IP states using lattice QCD. For the £,i(5830) and fi* 2 (5840) mesons, the results are in good 
agreement with the experimental values. Two further mesons are expected in the quantum channels J p = 0 + and 1 + near the BK 
and B*K thresholds. A combination of quark-antiquark and B 1 *' meson-Kaon interpolating fields are used to determine the mass of 
two QCD bound states below the B { *^K threshold, with the assumption that mixing with /f v ! V/ and isospin-violating decays to B'^n 
are negligible. We predict a J p — 0 + bound state B s o with mass mj M = 5.711(13)( 19) GeV. With further assumptions motivated 
theoretically by the heavy quark limit, a bound state with nig sl = 5.750(17)(19) GeV is predicted in the J p = 1 + channel. The 
results from our first principles calculation are compared to previous model-based estimates. 
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Over the years experiments have uncovered a number of 
mesons involving heavy quarks that do not seem to fit the sim¬ 
ple quark-antiquark picture suggested by quark models. Exam¬ 
ples of these include states in the charmonium and bottomo- 
nium spectrum 0]] as well as the charm-strange Z)* 0 (2317) and 
D sl ( 2460) 0 . The latter states are identified with the j = \ 
heavy-quark multiplet, where j is the total angular momen¬ 
tum of the light quark J3j] . These were predicted to be broad 
states above thresholds in potential models 0,0,00]. How¬ 
ever, the observed Z?* 0 (2317) and 7)^(2460) are narrow states 
below the DK or D* K thresholds 0 ], and it has been suggested 
that the thresholds play an important role in lowering the mass 
of the physical states [ 0 ]. In a recent lattice QCD simulation 
00 [llll these states are seen as QCD bound states below 
threshold with a mass in good agreement with experiment. 

In the B s meson spectrum only two positive parity states 
are known from experiment 11121113111441 . the B s i(5830) and 
£*,,(5840). The LHCb experiment should be able to see the 
remaining two states (0 + and 1 + ), which are expected to de¬ 
cay into s-wave states by emitting either a photon or a n° [15]. 
On the theory side there are a number of phenomenological 


model and EFT mass determinations 
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300 [01 . a determination using Unitarized EFT based 
on low energy constants extracted from lattice QCD simula¬ 
tions jfj. and some lattice QCD calculations in the static limit 
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a prediction 0 taking into account explicitly only quark- 
antiquark operators and extracting only the ground states in 
the system. This strategy can lead to inaccurate results in the 
vicinity of thresholds where meson-meson scattering can have 
a significant effect. None of the previous lattice simulations 
clearly establish the states in question as either QCD bound 
states below threshold or resonances above threshold. It is this 
gap which we aim to fill with the current publication. 

In this letter we present results for masses of the p -wave 
states of bottom-strange mesons with spin and parity quantum 
numbers J p - 0 + , l + ,2 + . For the heavy-quark doublet with 
j p — | masses determined using only quark-antiquark opera¬ 
tors agree with those of the observed £ v |(5830) and £* 2 (5840). 
This, as well as calculated mass differences between heavy- 
light mesons, verifies our computational setup. Then we simu¬ 
late B (t) K scattering in the scalar (axial) channel and extract the 
scattering matrix. Bound state poles are found below threshold 
and their location determines the masses of the B s q and B s \. 

The gauge configurations are from the PACS-CS collabo¬ 
ration li3lH . They have 2+1 flavors of dynamical quarks 
(up/down, strange); the bottom quark is implemented as 
a valence quark. The light and strange quarks are non- 
perturbatively improved Wilson fermions. The lattice spacing 
is 0.0907(13) fm and the Pion mass is 156(7)(2) MeV. The 
lattice size is 32 3 x 64 and we use stochastic distillation [32] 
for the quark propagation as in our analysis of the D s mesons 
1191 llOL 11]. This allows to include contributions with annihila¬ 
tion diagrams. Further details including the u, d, and s quark 
parameters can be found in 0. 

The dynamic strange quark mass and the associated hop- 
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Table 1: Parameter values in the dispersion relation (TJ for both the B and 
B* meson in lattice units. For our uncertainty estimates we also use alternate 
parametrizations. 


Table 2: Selected mass splittings (in MeV) of mesons involving bottom quarks 
compared to the values from the PDG 01. A bar denotes spin average. Errors 
are statistical and scale-setting only. 



B 

B * 

Mi 

1.5742(16) 

1.5960(27) 

m 2 

2.16(29) 

2.21(43) 

m 4 

1.4(2.6) 

1.05(77) 


ping parameter k s used in 13111 differs significantly from the 
physical value. We therefore use a partially quenched strange 
quark m v s al ± m s s ea . Different determinations agree very well 
and yield the value for k s [0 which leads to the Kaon mass 
m K = 504(1)(7) MeV. 

The bottom quark is treated as a valence quark and the Fer- 
milab method 133. 34ll is used. See Ref. 135. Tol l for details 


of our implementation. In the simplified form that we use 
136, 37], only the bottom quark hopping parameter k/, is tuned 
non-perturbatively, while the clover coefficients ce and cb are 
set to the tadpole improved value ce — cb = c® = 1/wjL where 
uq denotes the average link. There are several ways of setting 
Mo and we opt to use the Landau link on unsmeared gauge con¬ 
figurations. Within this simplified approach the static mass M\ 
may have large discretization effects but mass differences are 
expected to be close to physical {381] and can be compared to 
experiment. Determining the bottom quark hopping parameter 
translates into determining the spin-averaged kinetic mass M 2 
of IS B j mesons from the lattice dispersion relation 137] 


E(p) = M\ + 


2 M 2 


cr W 4 v-i 4 


2s2 


(P 2 ) 

8 Ml 


( 1 ) 


where p = yT| for a given spatial extent L. After trying 
multiple forms a simplified form without a W 4 term is taker]/] 
and for the value Kb — 0.096 used in our simulation we ob¬ 
tain M^y — 5086(135)(73) MeV. This value is significantly 
smaller than the physical value (ms, + 3m^)/4 = 5403.2/]g 
MeV but the effects on the binding energies used in our anal¬ 
ysis are small. This can be seen from the moderate difference 
between D s and B s binding energies we obtain and will 
be accounted for in the systematic uncertainty. For the analy¬ 
sis of the phase shifts the dispersion relations for the Kaon ( K ) 
and the heavy meson (B or B*) are needed. For the heavy B 
mesons we again take Eq. ([]} with W 4 — 0 and the results are 
tabulated in Table |T| For the Kaon the relativistic dispersion 


relation Ex(p) — + P 2 is used. 

The discrete energy levels for our combined basis of quark- 
antiquark and B (t) K operators are extracted from time correla¬ 
tions using the variational method 139l 40l 41,42]. For a given 


Mass splitting 

This work 

Experiment 

m B * 

- m B 

46.8(7.0)(0.7) 

45.78(35) 

m B s * 

- m Bs 

47.1(1.5)(0.7) 

48.7/“ 

m B s 

- m B 

81.5(4.1X1.2) 

87.35(23) 

lUy 

~m m 

44.2(0.3)(0.6) 

62.3(3.2) 

2 m B 

~ m hb 

1190(11)(17) 

1182.7(1.0) 

2m K 

-- ,n Tt 

1353(2)(19) 

1361.7(3.4) 

2 m Bc - 

m „ b - "V 

169.4(0.4)(2.4) 

167.3(4.9) 


quantum channel one measures the Euclidean cross-correlation 
matrix C,y(f) = {Oi(t)0\( 0)) between several operators living 
on the corresponding time slices. The generalized eigenvalue 
problem disentangles the eigenstates | n). From the exponential 
decay of the eigenvalues A n {t) ~ exp (-E n (t - to)) one deter¬ 
mines the energy values E„ of the eigenstates by exponential fits 
to the asymptotic behavior. The overlap factors give the 

composition of the eigenstates in terms of the lattice operators. 
In order to obtain the lowest energy eigenstates and energy lev¬ 
els reliably one needs a sufficiently large set of operators with 
the chosen quantum numbers. All error values come from a 
jack-knife analysis, where the error analysis for the phase shift 
includes also the input from the dispersion relation ([]}. 

To test our heavy quark approach we calculate a number 
of mass splittings involving heavy-light and/or heavy-heavy 
mesons, see Table [2] The quoted uncertainties are statistical 
and from scale-setting only and the values are not intended to 
be precision results. In particular our lighter than physical bot¬ 
tom quark mass strongly affects the spin-dependent splittings, 
but the effect tends to cancel with discretization errors. Esti¬ 
mates for both sources of uncertainty will be taken into account 
in our prediction of B s mesons. 

Partial wave unitarity implies that the scattering amplitude 
T(s) for elastic B { ' ] K scattering can be written as 

y[sT~ l (s) = p cot 5(s) - ip , (2) 


where p{s) is the momentum and s — E 2 the energy squared in 
the center of momentum system. Assuming a localized interac¬ 
tion region smaller than the spatial lattice extent a relation be¬ 
tween the energy spectrum of meson-meson correlators in finite 
volume and the infinite volume phase shift 6 has been derived 

| IJ El El [45, 


f(p) = P cot 6(p) = 


2 Zoo(i;(f) 2 ) 1 1 2 

-—p-- ~ — + ~F)P 

L yn ao 2 


(3) 


'The determination of the kinetic mass Mi (including its uncertainty) and 
thereby what is identified with the "physical” meson mass is rather insensitive 
(i.e. varies by < 15% of the uncertainty) to including or not including a VV 4 
term. (This is not the case for M\ and its uncertainty.) 


which applies in the elastic region and in the rest frame. Z .00 
denotes the generalized zeta function [44, 45] This real func¬ 
tion f(p) has no threshold singularity and the measured values 
can be found indeed above and below threshold. For ,v-wave 
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Figure 1: Plot of ap cot 8{p) vs. ( ap ) 2 for BK scattering in 5-wave. Circles 
are values from our simulation; red lines indicate the error band following the 
Liischer curves (broken lines). The full line gives the linear fit 0 to the points. 
Below threshold \p\ is added and the zero of the combination 0 indicates the 
bound state position in infinite volume. Displayed uncertainties are statistical 
only. 


scattering an effective range approximation (see Eq. (|3}) may 
be used to interpolate between the closest points near thresh¬ 
old. The imaginary contribution to T -1 becomes real below 
threshold (responsible for a cusp in Re T ). When the two con¬ 
tributions cancel, T 1 (see Eq. ©I develops a zero where 

mp B \) + \Pb\ = 0 . (4) 


That zero below threshold corresponds to a bound state pole of 
T in the upper Riemann sheet. 

For J p = 0 + we computed cross-correlations between four 
sb (in the form given in Table XIII of [Tol l ) and three BK (irre¬ 
ducible representation A|) operators: 

Of, = 0 PK — [S75M] (p = 0) [ujfb] ( p - 0) + {u — > d) , 

0 6 = 0 PK = [ sy,y 5 u\ (p = 0) [uy,ysb] (p = 0) + [u -> d) , 

0~j = 0 PK = Yj [ s T5«] (p) [uy 5 b\ (~p) + [u -» d) , (5) 

P= ±e x,y,z 2 lr /t- 


where we assume that the closeness of the BKn threshold can be 
ignored for our simulation. All operators are built according to 
the distillation method from quark sources that are eigenvectors 
of the spatial Laplacian, providing a smearing with a Gaussian- 
like envelope. The gauge links are four-dimensional normalized 
hypercubic (nHYP) smeared |47]. 

We omit B\ n interpolators since we work in the isospin limit 
where such decays cannot occur. We also neglect B[*\j, par¬ 
tially motivated by the threshold lying <9(140 MeV) above the 
B^* ] K threshold. Inclusion would necessitate a coupled chan¬ 
nel study which would need several volumes and considerably 
complicate the calculation. 

As in earlier experience it turned out that the full set of op¬ 
erators gave noisier signals than suitable subsets so for the final 
analysis we use the operator set (1,2,4,5,7). The energy values 
resulting from correlated 2-exponential fits to the eigenvalues 
are given in Table [3] 

In this channel B and K are in s-wave. If there is a bound state 
one expects an eigenstate with energy approaching the bound 
state energy from below in the infinite volume limit. The levels 


Table 4: Systematic uncertainties in the mass determination of the below- 
threshold states with quantum numbers J p = 0 + ,l + . The heavy-quark dis¬ 
cretization effects are quantified by calculating the Fermilab-method mass mis¬ 
matches and employing HQET power counting 0 With A = 700 MeV. The 
dominant contributions arise from mismatches in mg and niE and their size as 
a fraction of the reference scale A can be seen in Fig. 3 of 0. The finite 
volume uncertainties are estimated conservatively by the difference of the low¬ 
est energy level and the pole position (see also Equations (9) and (28) of [4811). 
The last line gives the effect of using only the two points near threshold for the 
effective range fit. The third point might be affected more strongly by the B { * ] q 
threshold, so it is reassuring that the difference in results between two-point and 
three-point fits is minimal. The total uncertainty has been obtained by adding 
the single contributions in quadrature. 


source of uncertainty 

expected size [MeV] 

heavy-quark discretization 

12 

finite volume effects 

8 

unphysical Kaon, isospin & EM 

11 

b-quark tuning 

3 

dispersion relation 

2 

spin-average (experiment) 

2 

scale uncertainty 

1 

3 pt vs. 2 pt linear fit 

2 

total 

19 


above threshold then would be dominated by BK operators with 
back-to-back momenta. This is exactly what is seen from the 
overlap ratios: The lowest level is dominated by operators 1,2 
and 4, level 2 by the B(0)K(0) operator 5 and level 3 by the 
B(1)K(-1) operator 7. 

As shown in (0 we can use the values of p cotd'(p) from 
Liischer’s relation to determine the effective range parametriza- 
tion near threshold. The energy eigenvalues give the points 
shown in Fig. Q] together with a linear fit. The value and slope 
at threshold can be related to the scattering length and effective 
range: 


o,q K = -0.85(10) fm, rg K = 0.03(15) fm . (6) 

Equation 0 gives the bound state position. From this the bind¬ 
ing energy is estimated to be m/j+raf-mj,, =64(13)(19) MeV; 
thus, using the physical threshold as input to minimize system¬ 
atic effects, we predict a bound state B s o with J p - 0 1 at a mass 
of 


m Bs0 =5.711(13)(19) GeV. (7) 


The first error is due to statistics and the effective range fit, and 
the second value is our estimate for the systematic error with 
the main contributions due to heavy quark discretization, un¬ 
physical Kaon mass, and finite volume effects. Details of this 
uncertainty estimate are provided in Table [4] 

For J p = 1 + we computed cross-correlations between eight 
sb (in the form given in Table XIII of ifTotl l and three B*K (irrep 
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Table 3: Energy levels for J p = 0 + (upper set), 1 + (middle set) and 2 + (lower set). A correlated 2-exponential fit is used and m = + 3 m B *) with 

m = 1.62897(43) in lattice units, to denotes the reference point in the generalized eigenvalue problem. Energy 2 in the middle set corresponds to the 5^(5830). The 
lower set shows the naive energy level for the J p = 2 + and corresponds to the £*-,(5840) using the same operator basis used in fioh for the Z)* 9 . 


level 

^0 

basis 

fit range 

d.o.f 

Ea 

E - m [GeV] 

(ap) 2 

ap cot(<5) 

/r[GeV 2 ] 

p cot(<5) [GeV] 

1 

2 

01,2,4,5,7 

4-16 

0.53 

1.7735 (44) 

0.315 (9) 

-0.0128(19) 

-0.106(10) 

-0.0606(88) 

-0.231(23) 

2 

2 

01,2,4,5,7 

4-16 

1.05 

1.8213 (29) 

0.419 (6) 

0.0066 (13) 

-0.116(18) 

0.0312(62) 

-0.252(40) 

3 

2 

01,2,4,5,7 

3-13 

1.35 

1.9139 (59) 

0.620 (13) 

0.0535 (35) 

-0.045(76) 

0.2532(165) 

-0.097(166) 

1 

2 

03,4,6,9,11 

4-14 

0.67 

1.7919(51) 

0.353 (11) 

-0.0141(22) 

-0.113(11) 

-0.067(11) 

-0.246(25) 

2 

2 

03,4,6,9,11 

3-14 

0.85 

1.8255 (42) 

0.428 (9) 

- 

- 

- 

- 

3 

2 

03,4,6,9,11 

3-14 

0.54 

1.8395 (45) 

0.457 (10) 

0.0050(24) 

-0.142(49) 

0.024(11) 

-0.308(106) 

4 

2 

03,4,6,9,11 

3-14 

1.19 

1.9406 (50) 

0.677 (11) 

0.0566(31) 

0.021(67) 

0.268(15) 

0.046(145) 

1 

2 

01,2 

4-14 

0.43 

1.8357(51) 

0.450(11) 

- 

- 

- 

- 


T+) operators: 

09 = °T,k = [ J r5«] (P = 0) [uy k b] (p - 0) + [u -> d) , 

Oio = Off = [sy,y 5 u\ (p = 0) [uy,y k b] (p = 0) + [u -> d) , 
0 n = Off = ^ [sy 5 u] {p) [uy k b] (-p) + {u -» d) . 

f^=±e xy z 2 Jt/L 

Comparing various subsets of operators the most stable set was 
(3,4,6,9,11), where four energy levels could be determined (Ta¬ 
ble 0 . 

Based on the overlaps, levels 3 and 4 are dominated by opera¬ 
tors 9 (B*(0)K(0)) and 11 (ZT(l)Zf(-l)), respectively. The low¬ 
est energy level (dominated by operators 3 and 4) agrees with a 
bound state interpretation. A linear fit to the points correspond¬ 
ing to energy levels 1, 3 and 4 gives the scattering parameters 

af K = -0.97(16) fm , rf K = 0.28(15) fm . ( 8 ) 

This indicates a B*K bound state B s \ with a binding energy of 
71(17)(19)MeV. Using again the physical threshold as input we 
obtain 


m Bsl =5.750(17)(19) GeV . (9) 

This state has not (yet) been observed in experiments. 

Notice that our determination assumes that the effect of s- 
wave - d -wave mixing is negligible on the scale of our uncer¬ 
tainty. This is motivated theoretically by the heavy quark limit 
ifitl (where such mixing is absent), which should be a good ap¬ 
proximation for bottom-strange mesons. 

Level 2 (dominated by operator 6 ) lies just below threshold. 
This is interpreted, as in the case of the D J i(2536) lHoh . to be 
the 7 = | state with J 1 ’ — 1 1 which does not couple to B K 
in ,v-wave in the heavy quark limit 01. The composition of the 
state with regard to the qq operators is fairly independent of 
whether the B*K operators are included or not. Assuming that 
the coupling to B*K in s-wave is indeed small, the “avoided 
level crossing” region is so narrow that this state may be treated 
as decoupled from the B*K scattering channel. Taking the mass 
difference with respect to the B s spin average and adding the 
physical value gives 

m B ' sl =5.831(9)(6) GeV, (10) 



Figure 2: Plot of ap cot 6(p) vs. (ap) 2 for B*K scattering in s-wave, as given 
by the levels 1, 3 and 4 in Table[3] see analogous caption of Fig. |7] 


where the errors are statistical and scale-setting only. In ex¬ 
periments one finds a resonance ZLi(5830) decaying 

dominantly into B* + K~ 10 MeV above threshold at 5.8287(4) 
GeV. The masses are in excellent agreement. 

The lowest energy level with J p - 2 + (irrep 7'j j correspond¬ 
ing to the B*,(5840) is extracted using just sb operators. The 
resulting mass is 


m Bsl = 5.853(H)(6) GeV, 


( 11 ) 


consistent with the observed value 0 ]. 

In summary we have analyzed the spectrum of positive parity 
B s mesonsf] and find two bound states below threshold, corre¬ 
sponding to the as-yet-unobserved B * 0 and /i vl IP states. Ta¬ 
ble [3 compares our first-principles lattice QCD calculation to 
previous results. Different variants of Unitarized ChPT along 
with phenomenological or lattice input (in particular dlHi) 
lead to mass predictions that are in good agreement with our 
calculation. Also, the model based on heavy-quark and chiral 
symmetry by Bardeen, Eichten and Hill 115] gives results that 
are remarkably close. 


2 The binding energies of the corresponding D s mesons were also reanalyzed 
with our updated procedure (basis, dispersion relation, etc.) and are fully com¬ 
patible with our old results hi and, within systematic uncertainties, with 
experiment. 
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Table 5: Comparison of masses from this work to results from various model 
based calculations; all masses in MeV. 


J p 

0 + 

1 + 

Covariant (U)ChPT [24] 

5726(28) 

5778(26) 

NLO UHMChPT (19J 

5696(20)(30) 

5742(20)(30) 

LO UChPT [17,18] 

5725(39) 

5778(7) 

LO*-SU(3) [16] 

5643 

5690 

HQET + ChPT [20] 

5706.6(1.2) 

5765.6(1.2) 

Bardeen, Eichten, Hill [15] 

5718(35) 

5765(35) 

rel. quark model [5] 

5804 

5842 

rel. quark model [22] 

5833 

5865 

rel. quark model f231 

5830 

5858 

HPQCD [30] 

5752(16)(5)(25) 

5806(15)(5)(25) 

this work 

5713(11)(19) 

5750(17)(19) 
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